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ABSTRACT 

The development of high-degree interpolation poly- 
nomials which use the values of the function and its sub- 
sequent derivatives is discussed. It is shown that if data 
of this type are available, high-accuracy interpolation is 
possible under the restrictive conditions of large step- 
sizes and few data values. 



CONTENTS 

Abstract ii 

INTRODUCTION 1 

THE GENERAL INTERPOLATING POLYNOMIAL 2 

DETERMINATION OF THE COEFFICIENTS 3 

SAMPLE DERIVATION 6 

FORMULAS FOR N' 0, 1. 2; n 0, 1, 2 8 

NUMERICAL EXAMPLES 12 

CONCLUSION 12 

ACKNOWLEDGMENT 14 

References 14 



GENERAL LAGRANGIAN INTERPOLATION FORMULAS 


by 

C. E. Velez 

Goddard Space Flight Center 


INTRODUCTION 

A computing system requiring an accurate numerical interpolation process may have, in addi- 
tion to the function, information concerning its derivatives at some set of discrete points. For 
example, consider a multistep numerical integration process designed to solve an equation of the 
form 

y F(x, y, y ' ) . 

The data concerning the function y(x) available at some point during the computation could then 
be written as the array 


y( x o) 


y ' 

( x o) 


y ( x o) 


y( x o • 

h ) 

y ' 

( x 0 

■ h) 

y ( x o 

• h) 

y( x o ‘ 

2h ) 

y ' 

( x 0 

* 2h) 

y ( x o 

4 21 1) 

y( x o 1 

kh) 

y ' 

( x o 

4 kh) 

y"( x o 

+ kh) 


where h is the integration stepsize, and k the number of ’’back points” required by the process. 

If the interval h of integration is modified, so that values are required for y, y' and y , at some 
points x< (x 0 , x 0 + kh) , then an approximation of these required values could be computed by 
applying an interpolation polynomial of the form 

n 

( a ) P(x) a,(x) f(x„ * ih) n :I k 

i - 0 
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to each column of the array independently (with appropriate coefficients), or by differentiating the 
polynomial approximation of y to obtain values for y' and y " . However neither of these methods 
uses all the information available about the function near the point x; for example, the approxi- 
mation polynomial for y would disregard the values y' (x 0 + ih) and y"(x 0 + ih). Therefore, 
consider using an interpolating polynomial of the form 


(b) P ( x ) " ^ 3i(X) f ( X ° + iH ^ * bi( - X) f ’( X ° + ih ) + c i ( x ) f" (x 0 + ih) 

i =0 


to obtain an approximation for y. Using such a polynomial would increase the accuracy con- 
comitant with the use of an increased amount of information about the function near the point x. 
Also an accurate interpolation may be possible even in a situation where the number of values of 
the function (k) is insufficient, or the stepsize (h) is too large for an accurate interpolation 
using a polynomial of the form (a). In the following paragraphs, the development of polynomials 
of the form (b) will be discussed for the case in which an arbitrary number of successive deriva- 
tives are available at a set of discrete points. In addition, the applicability of such polynomials 
will be demonstrated numerically. 


THE GENERAL INTERPOLATING POLYNOMIAL 

Let f be an N- times differentiable* function over the real line and assume the following 
data over the interval [x 0 , x k ] are given: 


H n > for 



0, 1, 2, * * * k 

0, 1, 2, • • * N ' 


where N' < N,f< n) r f^ n )(x i ), and x i e jx 0 , x k ] . The problem to be discussed then, is the formulation 
of polynomials of the form 

k N ' 

p (n) (x) -- h i" ) f i a> ' a) 

j = o i = o 


satisfying 


p(n) (x i ) = f . (n) , i = 0, 1, 2, • • * k , 


( 2 ) 


with remainder K n (x) - f (n) (x) - P (n) (x), where xe[x 0 , x k ] . 

*It is assumed that higher order derivatives remain bounded in the region of interest. 
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Some well-known examples of polynomials of this type are: 

(1) N' = o, the Lagrangian interpolation polynomial of degree k - n with its coefficients 
given by 


h<"> 


(X) = 4/">(x) = [(x - x i )/(x j - Xj)]} 


( 3 ) 


where 



(2) N' l , the Hermite interpolation polynomial of degree (2k + 1) - n with its coefficients 
given by 


A n 

h( n ) ( x \ - — — 

0-> <X) 


dx n 


|l - 2(x - Xj ) >•' Xj - x .j- i, 2 


(X) 


and 


( 4 ) 


(«> ■ [(-- »jv <»>] 


where ' (x) is the Lagrange coefficient given by Equation 3 for n o and 


k 

E 

i - 0 

. } 


(3) k o and n o, the truncated Taylor expansion about the point x 0 . 

In the following paragraphs, a general scheme for obtaining the coefficients of polynomials of 
the type given in Equation 1 will be discussed, and explicit formulas of the Equations 3 and 4 types 
will be given for x ' n, l, 2 and n - o, 1, 2. 

DETERMINATION OF THE COEFFICIENTS 

Since (N' + 1 )( k + 1) data points are given, it is seen that the coefficients of Equation 1 
could be uniquely chosen so that P (n) (x) is a polynomial of degree k(N # + 1) + (N' - n). Consider 
the case n 0 and let ^ - k(N' - l) + n'. In this case, such a polynomial could be constructed 
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by solving the following determinantal equation for P(x) (Reference 1, pp. 33-34): 


POO 

fi 

fk 

fi 




i 

i 

i 

i 

0 


2x n 


2x. 


0 0. . . . - 1) O - N # )x k v ~< N ' + 1 > 


= 0 - 


( 5 ) 


This technique, however, even for n' = 1, involves formidable computations. An alternate method 
which simplifies this computation considerably is suggested by Householder (Reference 2, pp. 
194-195) for N' = 1 . This technique, extended to the case of arbitrary N', is as follows: From 
Equations 1 and 2, 


f . = E Z>-.K)', ro 

j=0 i = 0 


( 6 ) 


This relation must hold, whatever the values of ft 0 ; in particular, we may have f. = 5 . 
and f (i > - 0 for i - 1 , 2, . . . N' and all j , so that Equation 6 reduces to 

i 

f p ; 1 ; h o,p ( x p) : 

but, since p is arbitrary, 

h 0,j ( x j) = 1 - j = 0, 1, • • • k . ( 

1 , while f . ~ 0 for j / m , and f j (i) r o for i = 1, 2, . . . N' 

f P z 0 = h o, m ( x p) ; 

h o,j ( x <0 = & j. q - (7) 

Moreover, if for some n, f <"> = l , while all other f (i) = o , then 
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Also, if for some m / p, f m = 
and all j , then 

i.e., 


BIllilil'HI III 


in 1 i i 1 1 mi iii 1 1 i mu 


f = 0 = h 


w 


i.e., 


h i.j ( Xq ) " 0 for i = 1, 2, . . . N' . 


Combining Equations 7 and 8 gives 


( 8 ) 


K i.j ( Xq ) _ S 0.i S i,q • 

Repeating the foregoing argument for the polynomials P Cn) (x), n 1 , 2, ... N', we obtain 

h i.i W (*0 = S i.n S i.q ' (9) 


Now consider the equation 

h i, j ( x ) = <Pi ,j O) (N +1) ( x ) 


(10) 


where the t. (x) are the Lagrange polynomials of degree k, and <P itJ (x) is a polynomial of de- 
gree n ' ; note that h. . (x) is then a polynomial of degree v and satisfies Equation 9 for q / j . 
Using the fact that a polynomial of degree v satisfying Equation 2 is unique, reduces the problem 
to determining ^^(x) so that Equation 9 is satisfied for q 3 j . This can be done by determining 
the N' derivatives of } evaluated at the point x - Xj using Equations 9 and 10, and then form- 
ing the polynomial 


0 itj (X) 


( x - x i) n ( x J 

zL n! 


(ii) 


We will then have the required polynomial 


P(x) 



( 12 ) 


with the coefficients h t j given by Equation 10. Polynomials p(")(x) for n^o can then be found 
by differentiating Equation 12: 


P (n) (X) 


£'<■> ■ L L 


h\"] (X) 


f ( i > 

J 
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The error term for P (x) can be derived using Equation 5 and the mean value theorem pre- 
cisely as indicated in Reference 1, pp. 22-24 for the case n = o. This technique yields 


R(x) 




(v + iy. 


(13) 


where £ e [x 0 , x k ] . As before, error terms for P (n) (x) can then be obtained by appropriate 
differentiation of Equation 13, noting that f is a function of x. Details concerning this differen- 
tiation can be found in Reference 3, pp. 66-67. 


SAMPLE DERIVATION 

As an illustration of the foregoing, a polynomial of the type 


k 2 

p(x) = L H h; ' i(x) fj(i) 

j =0 i-0 


will be determined. From Equation 9, we have the conditions 


h o.i 

( x 0 = s i, q h i.j 

(Xq) - 0 

h 2 , j ( Xq ) = 0 

h V 

i ( x< 0 = 0 h 'l. 

j (*0 = S i. q 

h ViH = 0 > 

h "o.j 

(xq) = 0 h"j 

j H = 0 

h *a.i ( Xq ) = S i.^ J 


from which to determine polynomials of degree 3k + 2 of the form 

h i tj O) = 4> iti O) 4 3 (x) , 

where . . (x) are quadratic polynomials. Differentiating Equation 15, we have 


i ( x ) 

□ 4>' . 

i ( x > 4. 3 ( x ) 

+ 3 <£>. . 

1 . J 

(x) 4. 3 (x) ( x ) 

,( x ) 

= *"i. 

j (x) 4 3 (x) 

+ 6 <*>' it 

j (x) l 3 (x) 4'. (x) 


+ ( x )[ 6 ^j OK,' 2 ( x ) + 3 4. (x) V(x>] 


(15) 


( 16 ) 
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From Equations 14, 15, and 16 we see that 

h o,iW = *o.i N 


1 


also 


that is, 


also 


that is, 


that is, 


h l.i( x i) = 

* 1.1 ( x i) 

= 0 

h 2,j ( X i) = 

*2.1 ( X i) 

= 0 ; 

h o. j ( X i) = 

*i.i N 

+ ( x i) = 0 ’ 

*0.1 ( x i) = 

-3*i ( x i) 


h i,i N = 

*i.i ( x i) 

= 1 

h i., ( x , 

*i.i ( x .) 

= 0 ; 

h o,j( x j) = 

*0,1 ( x i) 

- 124; 2 ( Xj ) + 34; (*.) 

*0.1 ( x ,) = 

i2t ; a ( x i 

) _ 3*1* ( x i) ’ 

h M( x i) r 

*i.j ( x i) 

+ 64; ( Xj ) - o , 

**• i ( X i) = 

- 6 *; ( x i) 


h ;.i(*i) = 

*2,j ( X l) 

= 1 . 


Hence, using Equation 11, we have 

*o.i < x ) z 1 ' 3 ( x ■ x i)' e i( x i) + " 2 ^ [ 12V j 2 ( x j) - ( x i 


*1.1 < x ) = ( x " x i) _ 3 ( x " x i) 2 l \ ( x i) - 

*2,1 O) 


( X ~ *i ) 2 
2 


and the required polynomial is given by 


P(x) 


K Z 

22 r. (x) f i 0) 


i = 0 


V (X) 


FORMULAS FOR N’ = 0, 1, 2; n = 0, 1, 2 

Case KLjL-Q 
For n' - 0) 


(17) 


k 

P(x) = Y h 0 J (x) f. (Lagrange), 
i = o 


(18a) 


where 

h o,j( x ) = l \ O) = v ‘ [( x " x i)/( x i " x i)] . 


and 


P'(x) = J2 h' , (X) f, 

j - o 


(18b) 


where 


^0, 


(x) - V. (x) - 't.(x) 2 ' - x ^ 


and 


p" ( x > = h °’> (x) { ‘ • 

1=0 


(18c) 


where 


h 


o. j 


(X) 


(x) 


fr 


l. (x)J 2 


LL 




8 



Case N f ~ 1 


For N # - i, 


k 

p (x) = ^ h o,j ( x ) f j + h i,j ( x ) f j (Hermite) , 


where 


h i,j ( x ) " ,j ( x ) V ( x ) . 


and 

0 O j (x> = i - 2 ( X - Xj ) -e; (x,) , 

( x ) ' ( x - *,) • 

where ' / ( x ) is given in Equation 18b with x x.. By differentiating, we obtain 


where 


k 

P ' (x) " 2^ h 0 .j O) f j + h i,j ( x ) f j 

J -0 


h . ,j ( x ) = y if j ( x ) ^ j 2 ( x > ’ 


(19a) 


(19b) 


and 


where 


>o.j ( x ) 



( x - x i) ( x ,) 


2 5 ' 


(X - X.) 



> 1 ., ( x > 


1 + 2(x 



*0 


P" (X) 


2 Z h; >" (x) { > 

j =0 


h l. J ( x ) 


f' , 


h i. j (*) = y i , j ( x ) V (*) • 


(19c) 
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y 0 ,j O) 


2 



2' 


X 


1 




<i) v i ( x i) 



Tj.i (x) - 


42' ^ 


x i) 




Case i<P=2 


For N' = 2, 


P(x) 


h °.i (x) f . + h i.i < x > f ; 

j =0 


h 2.j ( x ) { 'l . 


(20a) 


where 

h i. j (x) - (x) If (x) , 


and 


<*> 0 ., (x) - 1 + 6 


(x - Xj ) 2 j[e;( Xj )] 2 - I ^ (xjJ- 3(x - Xj ) l\ (x,) 


where i] (x.) is given by Equation 18c with 


■I.i (X) z (x - X,) - 3 (x - Xj) 2 V ( Xj ) 


<*2., (X) 



By differentiating, we obtain 


k 

P' (x) ' h o.j f j + h l.j (X) f i + (X) f i ’ 

j -0 


(20b) 


where 


h ! T J (x) > ifj (x) (x) 


and 


‘i.i (*> 


(X) 


(X- 


+ * 1.1 (X) 
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where <P iti (x) is given in Equation 19a, and 


and 


where 


*».}<*> = 

12 (x - x/) 

V ( x j) ~ T ( 

Xj) - 31 ’ X 

^i.j (x) = 

1 ~ 6 (x - 

x i) ( x i) 


* 2 ., ( X ) = 

( x ~ x j) ; 




k 

P " ( x ) r 2^ h o.j f i 4 h i,j < x ) f j + h 2 , j ( x ) f j 

j =0 


(20c) 


h i,j ( x ) r >i, j ( x ) y 3 O) 


and 


>’i., Cx) 


(x) i 6 c?; (x) > ' 


(x -■*,) 


3 0 , 


(x) 





where the 


(x) and <0| tj (x) are given in Equations 20a and 20b, and 


^0. j 

<*) ' 

« «; a ( x .) 


(X) 

- 6 ( Xj ) , 

*' 2.1 

(X) 

1 . 


If the data values are equally spaced over the interval [x 0 , x n J , letting h - - x._ 1 and 

s x - x p /h, the coefficients in Equations 18a, 19a, and 20a assume the following forms: 


(x) ‘ 

(“l) nj n ' ( s - i ) 
J ! (n - j ) ! 

(18a’) 

(X) = 

1 - 2 ( s - j ) S' j ± T ’ 

(19a*) 

(X) - 

h ( s - j ) , 
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4 > 0 , j ( x ) - l + 2 ( s - i Y 




+ S' 


0 - 0 s 


- 3 (s - j ) 2' 


(j " i) 


(20a>) 


4>i,l (x) - h[(s - j ) - 3 (s - j ) 2 2' (T^T)] 


h 2 , 

^2,j ( x ) " 2 ( s " J ) • 

NUMERICAL EXAMPLES 

In order to illustrate the power of the interpolation formulas including derivatives, Equations 
18a’, 19a’, and 20’ were applied to the following functions. 

(a) The solution of the initial value problem: 

y" = ■fijffjT - y(°) = y 0 ■ y'(°) = y ' 0 . 

where 

y = (y i - y 2 . y 3 ) 

and 

II y II H (y i 2 + y 2 2 + y 3 2 ) 1/2 • 


(b) The function y = sin (x) . 

In each case, the range of h was taken to be l/2 m where m = o, 4(1) and k = 1 , 10(3). The 
calculations were performed on the UNIVAC 1107 computer using double-precision arithmetic. 
These results are given in Table 1. 


CONCLUSION 

The development of interpolating polynomials utilizing any number of successive derivatives 
of the function has been discussed. It has been shown that for any fixed N' , such polynomials can 
be readily formulated. From the numerical results, we see that under the restrictive conditions 
of few data values and/or large stepsizes, accurate interpolation is possible by incorporating 
information concerning one or two derivatives. In general, it can be expected that for processes 
requiring accurate interpolation, if data concerning the derivatives are available or readily 
attainable, the polynomials discussed herein offer a distinct advantage of increased accuracy. 
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Table 1 


Numerical Results. 




Error st for example (a) 

Errors for example (b) 

m 

k 

Eq. 18a T 

Eq. 19a' 

Eq. 20a' 

Eq. 18a T 

Eq. 19a' 

Eq. 20a' 

0 

10 

.3 x 10 " 4 

.4 x 10' 8 

.2 x 10' 12 

.2 x 10 4 

.1 x 10' 13 

.4 x 10' I4 * 


7 

.2 x 10' 3 

.3 x 10' 6 

.8 x 10' 9 

.2 x 10 ' 3 

.2 x 10' 10 

.9 x 10' 15 * 


4 

.1 x 10' 2 

.3 x 10' 4 

.1 x 10' 6 

.5 x 10' 2 

.4 x 10' 6 

.8 x 10' 12 


1 

.7 x 10' 1 

.4 x 10' 2 

.1 x 10' 3 

.5 x 10 ' 1 

.1 x 10' 2 

.1 x 10 - 4 

1 

10 

.2 x 10 ' 5 

.1 x 10' 9 

.5 x 10' 13 

.3 x 10 ' 7 

.7 x 10' 15 * 

.9 x 10' 1S * 


7 

.5 x 10' 4 

.1 x 10' 7 

.4 x 10' 11 

.3 x 10' 7 

.3 x 10' 14 

.3 x 10' 14 


4 

.3 x 10' 3 

.3 x 10 6 

.6 x 10' 9 

.1 x 10' 3 

.4 x 10' 9 

.1 x 10' 14 


1 

.2 x 10' 1 

.2 x 10' 3 

.G x 10 s 

.7 x 10' 2 

.4 x 10' 4 

.8 x 10' 7 

2 

10 

.3 x 10' 6 

.5 x 10' !2 

.5 x 10' 14 * 

.lx 10 10 

.5 x 10' 14 

.7 x 10' 14 


7 

.2 x 10' 5 

.1 x 10' 10 

.1 x 10 ' 13 

.1 x 10' 7 

.1 x 10' 14 

.1 x 10' 14 


4 

.3 x 10' 4 

.4 x 10 7 

.2 x 10 11 

.9 x 10 ' s 

.2 x 10' 12 

.1 x 10' 14 


1 

.3 x 10' 2 

.2 x 10' 4 

.4 x 10 " 7 

.9 x 10' 3 

.1 x. 10 ' 5 

.G x 10' 9 

3 

10 

.1 x 10' 9 

.8 x 10' 14 * 

.1 x 10' 14 * 

.1 x 10' 13 

.1 x 10' 14 

.2 x 10' 14 


7 

.3 x 10 7 

.1 x 1 0 ' ' 3 

.2 x 10' ‘ 4 * 

.2x10' 10 

.7x10' ls * 

.1x10' 14 


4 

.2 x 10' 5 

.5 x 10' 11 

.4 x 10' 14 * 

.3 x 10' 6 

.1 x 10‘ 1S * 

.4 x 10' 1S * 


1 

.7 x 10' 3 

.7 x 10' 6 

.3 x 10' 9 

.1 x 10 3 

.4 x 10 7 

.5 x 10' 11 

4 

10 

.2 x 10' 11 

.1 x 10' 14 * 

.6 x 10' 15 * 

.4 x 10' 15 * 

.9 x 10' 15 

.1 x 10' 14 


1 7 

.2 x 10' 9 

.3 x 10' 14 * 

.2 x 10' 1S * 

.5 x 10' 13 

.4 x 10‘ 15 * 

.6 x 10' 15 * 


4 

.3 x 10' 7 

.2 x 10' 13 

.7 x 10' 15 * 

.1 x 10' 7 

.1 x 10' 1S * 

.2 x 10' 15 * 


| 1 

.4 x 10' 3 

.2 x 10' 7 

.5 x 10' 10 

.1 x 10' 4 

.1 x 10' 8 

.4 x 10' 13 

- 











* The error tabulated is the norm of the error vector. 

* Indicates the error within the accuracy of the analytic solution. 
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